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Binary particle coagulation can be modelled as the repeated random process of the combination of 
two particles to form a third. The kinetics can be represented by population rate equations based on 
a mean field assumption, according to which the rate of aggregation is taken to be proportional to the 
product of the mean populations of the two participants. This can be a poor approximation when the 
mean populations are small. However, using the Poisson representation it is possible to derive a set of 
rate equations that go beyond mean field theory, describing pseudo-populations that are continuous, 
noisy and complex, but where averaging over the noise and initial conditions gives the mean of the 
physical population. Such an approach is explored for the simple case of a size-independent rate 
of coagulation between particles. Analytical results are compared with numerical computations 
and with results derived by other means. In the numerical work we encounter instabilities that 
can be eliminated using a suitable 'gauge' transformation of the problem [P. D. Drummond, Eur. 
Phys. J. B38, 617 (2004)] which we show to be equivalent to the application of the Cameron- 
Martin-Girsanov formula describing a shift in a probability measure. The cost of such a procedure 
is to introduce additional statistical noise into the numerical results, but we identify an optimised 
gauge transformation where this difficulty is minimal for the main properties of interest. For more 
complicated systems, such an approach is likely to be computationally cheaper than Monte Carlo 
simulation. 



I. INTRODUCTION 

The process of coagulation or aggregation is responsi- 
ble for the coarsening of a size distribution of particles 
suspended in gaseous or liquid media. In essence the 
phenomenon consists of a sequence of statistically inde- 
pendent events where two (or possibly more) particles 
unite, perhaps as a result of collision, to create a com- 
posite particle, and each event reduces the number of 
particles in the systeirP^. This has further consequences 
such as colloidal precipitation or accelerated rainfall from 
cloudiPHSl Fragmentation can be present as welP, but 
precipitation processes are typically dominated by irre- 
versible agglomeration. The phenomenon is familiar and 
yet it can present some surprises, an example of which 
was presented by LushnikovS in an exact study of coag- 
ulation kinetics driven by various choices of aggregation 
rates. In the later stages of a process where the binary 
coagulation rate is proportional to the product of the 
masses of the participants, for example, a single parti- 
cle emerges with a mass representing a sizable fraction 
of that of the entire systerrP. This has been termed a 
gelation event and standard kinetic models of coagula- 
tion are unable to account for the phenomenon, for the 
simple reason that they are designed to describe systems 
consisting of very large populations of particles in each 
size class. They rely on a so-called mean field approxima- 
tion, though this is not always clearly recognised. When 
small particle populations play a key role in the kinet- 
ics, different approaches become necessary^, the most 
common being Monte Carlo modelling^. 

This paper investigates the utility of a rate equation 
model of kinetic processes that is able to capture small 
population effects. The Markovian master equations that 



describe coagulation may be transformed mathematically 
into a problem in the dynamics of continuous stochastic 
variables acted upon by complex noise. The solutions to 
the transformed dynamics may be related to the evolving 
statistical properties of the populations in the physical 
system. The recasting of the problem into one that con- 
cerns complex 'pseudo-populations' can be done in two 
distinct ways, either using methods of operator algebra 
similar to those employed in quantum field theorie a 12 " 22 ', 
or through the so-called Poisson representation of the 
populations^^. The physical problem concerning the 
stochastic evolution of a set of (integer-valued) popula- 
tions is replaced by the task of solvi ng an d averaging a 
set of stochastic differential equation d 27 ' 28 '. 

The advantage of such a tra nsform ation over Monte 
Carlo simulation of the processors' emerges for cases 
where the particles come in many species or sizes, since 
the configuration space of the populations, and hence 
the computational cost increases very rapidly, but in or- 
der to illustrate the approach we study a very simple 
example of coagulation, where the aggregation rates do 
not depend on the masses of the participants. This is 
very different from the cases that exhibit gelation, stud- 
ied by Lushnikov and others, and we might not expect 
major deviations from a mean field approach, but nev- 
ertheless it is possible to demonstrate the analytic pro- 
cedure, and compare the accuracy of numerical pseudo- 
population dynamics and averaging with respect to other 
approachesS2tt25| in order to form a view on the usefulness 
of the approach. 

In section [II] the master equations describing the basic 
problem of A + A — > A aggregation are developed and 
then the Poisson representation is introduced and used 
to derive the equivalent stochastic dynamics problem and 
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associated averaging scheme. The exact solution for the 
evolving complex pseudo-population is obtained and its 
properties established. In section |III| a parallel numeri- 
cal study of the stochastic problem is discussed. Inher- 
ent instabilities in the dynamics may be eliminated using 
the Cameron-Martin-Girsanov formula, or equivalently 
through a so-called gauge transformation introduced by 
Drummond 25 . However, this comes with the introduction 
of diffusive dynamics for a subsidiary variable that for ex- 
tended evolution times would seem to require high com- 
putational cost for accuracy. Nevertheless, with a certain 
choice of transformation, termed an optimised gauge, we 
can ensure that the coagulation process is completed be- 
fore this diffusive noise becomes apparent. We comment 
on the procedures and prospects for their further use in 
more complicated models in section |TV} 



and this can be integrated to give 



II. ANALYTIC MODEL OF COAGULATION 

A. Master equations 

We consider the dynamics of a population of particles 
of a species A undergoing binary coagulation A + A — > 
A. The distinction between particles of different mass is 
ignored and the rate at which aggregation events takes 
place is a constant, making this one of the simplest cases 
to study. The evolution of P(N,t), the probability at 
time t that N particles remain, is described by a set of 
master equations of the form 



dP(N,t) 
dt 



= k(N + 1)NP(N + 1, t) - kN(N - 1)P(N, t). 

The first term corresponds to the rate of gain in the prob- 
ability of observing population TV as a result of an ag- 
gregation event in circumstances where there are N + 1 
particles, properly weighted by the number of particle 
pairs in such an initial state and the probability re/2 of 
an event per unit time and per pair. The second term 
corresponds to the rate of loss of probability due to ag- 
gregation starting from a population equal to N. Multi- 
plying the master equations by N and summing over N 
it can be shown that 



dN 
~dt 



N 2 - N 



(2) 



where statistical averages are written F(N) — 
£~ F(N)P(N,t). _ 

If we take the view that TV 2 N 2 , which corresponds 
to the neglect of fluctuations in population (which is why 
the procedure is called a mean field approximation), and 
consider a large mean population N 3> 1, such that only 
the first term on the right hand side in Eq. ^ is retained, 
we can write 



dN 
~dt 



-kN 2 , 



(3) 



N 



Nn 



1 + nN t 



(4) 



where N is the initial mean population. This is the clas- 
sical Smoluchowski solution to this type of coagulation 
kinetics^. 

If the second term in Eq. ^ is retained, but again 
neglecting fluctuations, then the integration yields 



N 



No exp(ret) 



1 - N (1 - exp(ret)) ' 



(5) 



which ought to be more accurate than Eq. Q , especially 
in the limit t — > oo when it goes to unity rather than zero. 

Rather than making the mean field approximation 
TV 2 « N , we could generate an evolution equation for 
N 2 by multiplying the master equation by N 2 and sum- 
ming, namely 



dN 2 
dt 



«(■ 



■27V 3 + 37V 2 - N 



(6) 



but now the third moment appears on the right hand 
side. An equation for the evolution of this moment would 
involve the fourth moment: this is a hierarchy problem 
that often arises in kinetic theory. A closure condition 
such as N 3 = N 2 N could be imposed, but the accuracy 
of such assumptions is questionable. 

Similarly, the numerous master equations for the 
P(Ni, Ns, . . . , t) in a more general model where there 
are different categories of species {At} with populations 
{Ni} would reduce under a mean field approximation to 



dNj 
dt 



1 i_1 



-3,3 



N, 



(7) 



where the coefficients re, 



xj.j quantify the rate of aggrega- 
tion of the form Ai + Aj — > Ai + j between species i and 
j, with i representing species mass, for example. If we 
should need a treatment beyond the mean field approxi- 
mation, then for relatively simple stochastic processes we 
could use Monte Carlo simulations to extract the relevant 
statistical properties, but for more complex coagulation 
problems the use of this technique can become computa- 
tionally expensive. 

The transformation of the problem we wish to investi- 
gate is of interest since it reduces the master equations to 
the form of Eq. Q and more generally Eq. ^ , but sup- 
plemented by a noise term on the right hand side. This 
is not the same as inserting noise to represent a random 
source term in the population dynamics, nor is it an ad- 
ditional representation of stochasticity in the coagulation 
kinetics. It is a noise term that creates the statistical cor- 
relations in the populations that are neglected when the 
mean field approximation is taken. It turns out that the 
dynamics of simple binary coagulation may be modelled 
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by a pseudo-population 
equation 

d(f> _ 
~dt 



t hat evolves according to the 



(8) 



where £ is a noise with certain statistical properties. This 
may be compared with Eq. ([3]). The notable aspect of 
this representation is that the noise is complex, such that 
the variable <j> is generally complex as well. Averages are 
to be taken over the noise history to make a connection 
with the dynamics of a real population. Nevertheless, the 
solution of Eq. Q can be a simpler task than that posed 
by Eq. (|T|). We next describe how this transformation 
may be achieved. 



B. Poisson representations 

Gardiner and ChaturvedP^ outlined a method for 
transforming master equations into a Fokker-Planck 
equation by representing probability distributions as in- 
tegrals of weighted complex Poisson distributions. The 
probability of finding N particles at time t may be writ- 
ten 



P(N,t) = (h f{4>, t)exp(- 



(9) 



This is a superposition of Poisson distributions, over a 
closed contour of complex mean values <f>, with an evolv- 
ing weighting function f((f>, t) which has initial value 

/(0,O)=^P(7V,O)^0- Ar - 1 exp(0), (10) 



Substituting the Poisson representation of the proba- 
bilities given in Eqs. (|9| or (13 1 into Eq. ([!} leads to an 
evolution equation for f((j>,t): 



0/ d(£f)_ d 2 (0 2 f) 
dt K d<j> K 84> 2 ' 



(14) 



This emerges as long as integration boundary terms can 
be dropped, which could potentially be a problem for the 



2-d integration scheme ( 13 1 but is not an issue for the 
closed contour integral representation |9p31 But since 
a stochastic problem described by a Fokker-Planck equa- 
tion such as ( [Til ) can b e recast as one involving a stochas- 
tic differential equation (SDE) , the properties of the dis- 
tribution f(4>,t) can be reconstructed by solving 



-K(j) 2 dt + i(2 K ) 1/2 (j)dWt, 



(15) 



where dWt is an increment in a Wiener process with 
properties (dWt) = and ((dWt) 2 ) = dt, the brackets 
representing an average over the noise. This Ito-type 
SDE takes the promised form of Eq. Note that 

the noise term is complex since the second derivative in 
Eq. (14 1 has a negative coefficient. The equivalence be- 



tween approaches is such that the average of a function 
A(<j>) weighted according to the solution f(4>,t) with ini- 
tial condition f(<f>, 0) = 5((f> — 4>o) (namely the Green's 
function of the Fokker-Planck equation) is equal to the 
average of the same function of a stochastic variable cj>(t) 
evolved over a noise history Wt with initial condition 
(f)(0) = 0o- Both averages will be denoted with angled 
brackets and subscript in the form (A((j)))^ with time 
dependence understood. We next demonstrate that <p(t) 
can be determined explicitly for this problem. 



such that if the contour C includes the origin, the initial 
condition P(N, 0) is recovered from Eq. (|9|. There are 
two particular initial conditions of interest. If the initial 
population is known to be Nq then P(N,0) = 5nn and 

f(<j>,0) = f N M = ^0- Ar »- 1 ex P (0), (11) 

whereas if P(N, 0) is a Poisson distribution with mean A 
then we might use 

/(&o) = /f(0)=£^i£r (12) 



2ni <j) N + 1 



Alternatively, we could employ a Poisson representation 
that involves a 2-d integration over the entire complex 
plane, namely 



lN 



P(N,t)= I f(4,,t)exp(-4>)^d^ 



(13) 



which is particularly convenient if P(N, t) is initially a 
Poisson distribution since we would then use f(<p, 0) = 



C. Stochastic evolution of a complex 
pseudo-population 

Consider G(<j>,x) — _1 exp(x) with stochastic vari- 
ables evolving according to dip = adt + bdW t and dx 
I nil + qdWt ■ By Ito's lemma we have 



dG 



so 



9G„ dG , 1,J 2 G, , d 2 G 

gJ^ + -d X+ -t 2 w dt +bq 



dt- 



1 o 2 G 



dt. 



d<j)dx 2 dx 2 
(16) 



I b 2 ba a 2 

dG = — Gd(f> + Gdx + —Gdt - —Gdt + —Gdt (17) 



q GdW t 



a b 2 bq q 2 \ 



With the aim of reducing this to a deterministic form, 
we choose b = qd>, in which case 



dG = {-4>^a +p + q 2 /2) Gdt. 



(18) 



We set q = i(2k) 1 / 2 and p = k, in which case the SDE for 
x integrates to give x(t) = K,t+i(2n) 1 / 2 Wt, having chosen 



initial condition x(0) = 0. Furthermore if a = —K(f> 2 then 
d(j) = -n(j) 2 dt + i(2K) 1/2 0dW t , (19) 
as desired, and 

dG = n4>Gdt = ne x dt = K e Kt+l{2K]U2Wt dt, (20) 
and this integrates to 

G(t) = G(0) + k J exp (ks + i(2n) 1/2 W^ ds. (21) 

SmceG = (l>- 1 (t)exp(Kt + i(2K) 1 / 2 Wt) and G(0) = 

we obtain an explicit solution to the stochastic dynamics: 

<poexp(Kt + i(2n)iW t ) 



1 + K(f>o Jq exp(Ks + i{2k) 2 W s )ds 



(22) 



We next show that an average of this stochastic variable 
can be related to N, the particle population averaged 
over the stochasticity of the physical coagulation process. 



D. Integrals over initial pseudo-population 



Using Eq. ^ we can write 

AjN) = ^ A(N)j)f fat) exj>(- 



N=0 



1> N 

Nl' 



f{ct>,t)A{4>)dct> = {A{cj>)), 



(23) 



noting that the final expression differs from the average 
(A(0))0 O introduced earlier, since the initial weighting 
/(</>, 0) takes the general form of Eq. (10 1 instead of 
8((j} — 4>q), and where 



N 



Nl' 



(24) 



For example, if A(N) = N, then A(<f) — <fi, and for 
A(N) = N 2 we have A(<f>) = <t> 2 + <j). A similar con- 
struction might be made using the 2-d integration scheme 
( [13] ). As for the state probabilities P{N,t), we note that 
P[N',t) = E N SNN>P(N,t) = W-^ Taking A(N) = 
$NN' we find that the corresponding A((j>) according to 
Eq. JiiJ is exp(-(j))4> N ' /N'\ and hence 



P(N,t) = (eM-^)r/Nl). 



(25) 



As noted earlier, the Fokker-Planck equation (14 1, 
written as £/(</>, t) = 0, has a Green's function satis- 
fying CG(4>, <f>0i t) = for t > with initial condition 
G(</>, </>o, 0) = 5{(f> — 4>o) and we can write 



Im (j) ' 




Contour of initial cj)o 





Re0 



Figure 1. Sketch of stochastic realisations of the evolution 
of pseudo-population cj>(t) in the complex plane according to 
Eq. \15) , starting from point cj>o on a contour that encircles 
the origin. An average over all end-points, three examples of 
which are shown as filled circles, corresponds to the quantity 



where the contour should pass through the point <j) — 
4>o such that the initial condition (A(4>)),p = A(<f)o) is 
satisfied at t = 0. Furthermore, the evolution of f(cj>,t), 
starting from an initial distribution /(</>, 0) that is non- 
zero only at points cj>o that define an origin-encircling 
contour G, can be constructed from the superposition 



ffat) = J>dtt>oGfa<h,t)f(<h,0), 



(27) 

so that from Eq. fl23) : 

A(N) = (£,({, # O G(0, fo, t)f(4, , 0)i(<^»)#, (28) 



c Jc 



and we deduce using ( 26 ) that 

AjN) = $ # o /(0o,O)(i(<^ o , 



(29) 



which is intuitively understood as a superposition of 
the stochastic evolution of a function of the pseudo- 
population <f>(t) that evolves from points on a contour of 
initial values <fio, weighted by the function f(<f>o, 0), given 



in Eq. (10), and then averaged over the noise present in 
its SDE Eq. (IT5|) or in its solution (p2l. This picture is 



illustrated in Figure [T| 

Of course , we could equally well employ the rep- 
resentation A(N) = J f((j),t)A{(f))d 2 4> together with 
the associated 2-d integral versions of Eqs. (26 1 
and (27) to obtain an analogous result A(N) — 



j d 2 (f>of(4>Q,0)(A(<j))) ( f >0 which. is particularly convenient 
for a initial Poisson distribution of population, such that 
A(N) = (A((f>))\, where A is the initial mean. 



E. Noise-averaged pseudo-population 



(i(0)) 0o = £ G(0, 0o, (26) 



We now establish the statistical properties of the 
stochastic variable cj)(t). Our strategy will be to represent 
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{4>)cj, , the average of Eq. (22 1 over the noise for a given 
initial value (fio, as a (formal) series in positive powers 
of <f)Q. The first term is straightforward to identity: we 
write 



W*o« J dWt4> exp(Kt + i(2K)2W t ), (30) 
and then employ the stochastic integral identity 

dW t exp(icW t ) = (exp(icW t ) = exp(-c 2 t/2), (31) 
such that 

(<!>)*> =<h + 0{tft). (32) 



The proof of Eq. (31 1 follows from considering the 
stochastic evolution of F(x) — exp(x) where x = icWt- 
Ito's lemma states that dF = F'dx + (ic) 2 F"dt/2 = 
icFdWt - c 2 Fdt/2, where the primes indicate differenti- 
ation. Averaging leads to d(F) = —c 2 (F)dt/2 and hence 
to (F) = exp(-c 2 i/2) since F(0) = 1. 
The formal series is written 

oo oo 

<<«0o = E C i ®4 = <t>0 WWn(t) (33) 

j=l n=0 

with C 3 -(*) = (-K) J ' _1 M,-_ 1 (i) and where 



Kt+i(2 K ) 1/2 W t 



(34) 

We focus our attention on the evolution of the mean 
particle population N for cases where there are No parti- 
cles at t = 0. The initial distribution is P(N,0) = 5nn , 



so we use Eqs. (|29| and (J 1 1 [) to obtain 

" :i -1 exp((/> ) (</>)0 Q 



"00 7^7^0 



27Ti 



(35) 



and by inserting ( 33 1 and expanding the integrand we 
find that 



N 



N 



(36) 



using standard calculus of residues. Notice that the series 
for N is finite even though the series for (4>}^, is infinite. 
We have already established that C\(t) = 1 in Eq. 



(32 1. In order to evaluate C 2 (t) we need to consider 

I(t, s) = (exp (i^k) 1 ' 2 (W t + W.j) ) , (37) 
which in more explicit form reads 

r°° / 
I(t, s) = / dW t dW s eiq> U(2k) 1/2 (W t + W,)" 



1 



(2ir(t - s)) 1 / 2 
1 



exp 



(W t - Wsf 
2(t - s) 



(38) 



(27TS) 1 /: 



exp 



wl 

2s 



in which the weighting is a product of the gaussian prob- 
abilities for generating a value W s of the Wiener pro- 
cess at time s, and a value W t at the later time t given 
the earlier value. Writing Wt + W a — W[ + 2W S where 
W' t = Wt — W s we find this factorises as follows: 



I(t,s) = J dW^exp (i(2n) 1/2 W{ 



exp 



(2tt(<-s)) 1 /2 
x / dW s exp(2i(2K) 1/2 W ! 



W' 2 
2(t - s) 
1 



(2tt,s)i/2 exp 
(exp(i(2K) 1 / 2 W t - s )(exp(2i(2^ 1 / 2 W s ) 

e -n(t-s) e -ins _ e -K(t+3s) 



wl 

2s 



(39) 



using Eq. (31), and hence C 2 (t) — —K,M\(t) with 
Mi(t) = e Kt f dse KS (exp (z(2k) 1 / 2 (W t + W a ))) 



dse KS I(t, s) 



ds e 



-2ks 



2k 



(1 



- 2 Kt\ 



The averaging rapidly becomes more laborious. Consider 
next C^i) = K 2 M 2 (t). We need to evaluate 



M 2 = ( dsi / ds 2 e K(t+Sl+S2) 
./o 

^(2 K ) 1 / 2 (l¥ t + ^ Sl +VF S2 ) 

2 f dsi f 1 ds 2 e K( * +;il+S2) 
./o 



x ^exp 
ft 



'0 JO 

x (exp [z(2k) 1 / 2 (W t - W S1 + 2(W S1 - W S2 ) + 3W S2 )] ) , 

(41) 

where an ordering t > si > s 2 > has been imposed 
by the choice of integration limits, the change in which 
is accounted for by inserting the prefactor of two. Using 
Eq. (31 1 this reduces to 



M 2 = 2 f ds-L f ds 2 exp(K(i + s 1 + s 2 )) 
Jo Jo 

x exp (— n(t — si)) exp (— 4k(si — s 2 )) exp (— 9ks 2 ) 



2 / dsi / ds 2 exp(— 2ksi) exp(— 4ks 2 ), 
Jo Jo 



(42) 



which becomes 



M 2 



1 /•* 

— / dsi exp(— 2ks\) (1 — exp(— Ansi j) 
'k Jo 



2 
1 

2k 



2 k 



1-e 



- 2 k* 



-6nt\ 



(43) 



Similarly, we can show that 



M 3 = 6 f d Sl f 1 ds 2 f 2 ds 3 e- 2KSl e- 4KS2 e- 6KS:i , 
Jo Jo Jo 



(44) 
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and the pattern that emerges is 



M n = nl I dst 
Jo 



ds n e~ 2KSl ■ ■ ■ e - 2nKS " 



(45) 



In order to proceed further, we notice that the M n are re- 
lated to one another by repeated integration of the nested 
integrals. We define 

m n (Kt) = ^-M n {t) 



n< 

2nt psi fSn — 1 

ds\ I ds2 ■ ■ ■ I ds n e~ Sl e 
Jo Jo 



such that mi = 1 — exp(— 2nt) and 



Si „ — 2»2 



(46) 



is a difference of triangular numbers, such that 

H 



1 1 

n n[n + (n — 1)\ 



i 



-i) m [m^n+i 



+ (n|(2n+l-i) 



m (n), (51) 



where we define 

mo(n) = 1 — exp[— n(n + 
Next we note that 



(52) 



m 2 = \ (2 - 3e- 2Kt + e" 6 **) = I mi - 1 (l - e -<*««) . 
6 '26 

(47) 

We evaluate 7713 in detail in order to relate it to toi and 
to 2 : 



Jo Jo Jo 



/ dsi / ds 2 e- Sl e- 2s2 - (l-e- 352 ) 
./o 3 



1 1 f Sl 

= o TO 2 - n / ds i / ds 2 e" Sl e 
•J Jo 

1 1 f 2Kt 1 

= 3 7712 - 3 y ^-V 1 - 6 " 

1 11 111, _ 12 

= -ma mi H 1 — e 

3 3 5 3 5 6 V 



O 8 1 



(48) 



On the basis of an analysis of further cases, we conjecture 
that the pattern is 



Ij|(2n + l-j) 



3=1 



(2n + l-J)] = 



l\{2n)\ 
2 l {2n-l)V 



(53) 



so finally 



2'(2n -Q!(-l)' +1 
(2n)W. ' 



(54) 



in agreement with Eqs. (47 1 and (48 1, and where the 



n— dependence of too appearing in the final term is un- 
derstood. This provides a complete recursive solution to 
the stochast ic p roblem. Evaluation of the m n by integra- 
tion of Eq. ( 46 I using Mathematical for n up to twenty 



confirms this relation and the conjectured pattern that 
relates them. 



F. Mean populations in example cases 



1 m n -2 

m n = -run-! — — — r 

n n[n + [n — 1)\ 

™n-3 

n[n + (n - 1)] [n + (n - 1) + (n - 2)] 

mi 



+ (-1) 



n+l 



n [n + {n - 1)] • • • [n + (n - 1) H h 2] 

1 

n [n + (n - 1)] • • • [n + (n - 1) -I h 2] 



r-lnt 

x / c?si exp[— (n + [n — 1] 
Jo 



l) Sl }. 



(49) 

This may be simplified since n + (n — 1) + ■ ■ ■ + 1 = 
^n(n + 1), a triangular number, and 

n+ (n- 1) H h (n-i + 1) 

= n + (n - 1) H hi 

- [(n - 3) + (n - i - 1) + • • • + 1] 

11 j 
= -n(n + 1) - ~(n - j)(n - j + 1) = - (2n + l - j) , 

(50) 



From the previous section we can establish that 
Ci(t) = 1 and 

C 2 (t) = -\ mi = -\{l-e-^) 



(55) 



C 3 (t) = im 2 = i (2 - 3e- 2Kt + e" 6 **) 
^) = !-3 = ~^e--ie- + i Le-- 



and we use these and Eq. ( 36 ) to obtain the exact time 



dependence of the mean population for initial population 
No from 1 to 4. For No = 1, we obtain 



— No 1 



(56) 



which is clearly what we would expect. For Nq — 2 we 
find that 

3=1 ^ ° 



2-(l 



-2nt\ 



l + e 



-2Kt 



(57) 
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and this can be checked by solving the underlying mas- 
ter equations for this simple case. They are dP(l, t)/dt = 
2rcP(2,i) and dP{2,t)/dt = -2/cP(2,t) with initial con- 
ditions P(1,0) = and P(2, 0) = 1, giving solutions 
P(2,t)_= exp(-2Kt) and P(l,t) = 1 - exp(-2Ki) such 
that iV = £jv iVP(iV, t) = P(l,t) + 2P(2,<) = 1 + 
exp (— 2Kb). 

For iVp = 3 we get 



JV ! 



(iVo-j)! 



3Ci(i) + 6C 2 (t) + 6C 3 (t) 



j_ + 3 -2nt + 1 -6«t 



(58) 



and using Eqs. (55 1 we can construct the time deriva- 
tives, giving dCi(t)/dt = together with 



-2nt 



1 dC 2 {t) = 
k dt 

1 dC 3 (t) = 1 2Kt 

k dt 2 ^ 

1 dC A {t) = 3 c _ 2Kt 

k 20 



-6«t \ 



(63) 



-e 
4 



-12Kt 
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which may also be checked by solving the master equa- 
tions directly. 

Finally for Nq = 4 we obtain 



N ! 



3 = 1 



(N -j)\ 



4C 1 (t) + !2C 2 (t) + 24C 3 (i) + 24C 4 (i) 
1+ 9 e - 2Kt + e- 6K *+ 1 e" 12K *. 



(59) 



All these solutions satisfy the initial condition N = Nq 
at t = and tend to unity as b — > oo, as required. 

These results should be compared with those of 
Barzykin and TachiyaP^ who found an exact solution to 
the master equations for this problem using a generating 
function approach. They obtained 



Wo -1)1 iVp! .- J(J - 1)lrt ) 



V 1 'i-v,. ■ ./ i,:,v, 



, (60) 



and our calculations are consistent with this expression. 
We now go on to study the variance of the particle pop- 
ulation. 



G. Variance in population 



which puts us in a position to calculate the variance a 1 



Considering first the trivial case No — 1. we de- 
duce N 2 = 1 and a" = as we would expect since 
there is no evolution. For No — 2, TV 2 — 2 + 
2 [exp {-2Kb) - |(1 - exp {-2Kb))] = 1 + 3 exp {-2Kb). 

This has the correct limits of N 2 = 4 and 1 for b = 
and t — > oo, respectively. The variance is then a 2 = 
1 + 3exp(-2Ki) - (1 +exp(-2 K f)) 2 = exp(-2«;i) - 
exp (— 4Kb) . This also has the correct limits of zero at 
both t = and t — > oo. It goes through a maximum 
at —2k + 4«exp(— 2Kb) = or 2nt — In 2. The maxi- 
mum variance is then = 1/4. We can check the be- 
haviour of the variance for this case using the previously 
derived solutions to the master equation for No = 2. 
We have N 2 = J^n N 2 P{N,t) = P(l,t) + 4P{2,t) = 
1 — exp(— 2nt) + 4exp(— 2 Kb) = 1 + 3cxp(— 2Kb) as re- 
quired. 

We employ MathematicJ^ to calculate higher order 
coefficients m n and hence the C n {t) in Eq. (36 1. We il- 



lustrate this with the time-dependent mean and standard 
deviation for the case with No = 12, shown in Figure [2] 
The exact form of TV obtained from the analysis is 



In order to calculate th e va riance in the population we 
need the second moment N 2 , which is equivalent to the 
quantity (cf) 2 + <j)). We could expand <fi 2 as a formal power 
series, just as we did for (j>, but instead we exploit a short 
cut due to the relationship 



d{4>) 
dt 



(61) 



(361 



that arises from Eq. (15). We hence obtain from Eq. 



N 2 



N 

E 



k dt 



C 3 {b) 



N \ 



{N -j)V 



(62) 



N 



58786 



-132k* 



21e 



-llOrei 



209e 



-90«t 



+ 1309e 



-72Kt 



5775e 



-56k* 



19019e 



-42ret 



48279e 



-30/ci 



95931e 



-20/ci 



149226e~ 



177650e" 



149226e~ 



58786 



(64) 

which is consistent with the Barzykin and Tachiya 
expression^. Eq. Q happens to account reasonably 
well for the evolution of the mean population, as shown 
in Figure [2j but of course the mean field approximation 
upon which it is based cannot produce the standard de- 
viation, which is also shown. 
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N 

Eq. (5) 




Figure 2. Mean (solid line) and standard deviation (short 
dashed line) of the population as a function of time, in a co- 
agulating system with initial population No = 12 and with 
k = 1. The long dashed curve is the mean population accord- 
ing to Eq. (5}. 



H. Averaged pseudo-population as t — > oo 

It is of interest to determine the value of (4>}^> as t 
oo. Consider, using Ito's lemma, 



= n^^dcf) + ^(-2n<j) 2 )n{n - l)<j> n - 2 dt, 



and insert Eq. ( 15 ) and average such that 

d(<p n ) = -Kn((/) n+1 )dt - nn{n - l)(4> n )dt. 



(65) 



(66) 



Assuming that as t — >• oo all the (4> n ) become time inde- 
pendent, we deduce that (<fi n+1 ) = — (n — 1 )(</>") in this 
limit, in which case all moments for n > 2 tend to zero, 
for all initial conditions. 
Next, we note that 

d(e" ) = -e-^dcf) - K<f) 2 e-^dt = -i{2K)<f>e-^dW u (67) 



using Eq. (15 1, such that d(exp(— <fi)} = 0, and imposing 



the initial condition we obtain (exp(— (j)))^ = exp(— 4>o) 
for all t. Since 



E 

n=0 



(-1)" 



as t — ^ oo , due to the vanishing of moments with n > 2. 
we conclude that the stochastic dynamics of tp(t) give rise 
to 



(<P)4> -> 1 - exp(-0 o ) 



(69) 



in the late time limit . We shall use this result to check the 
outcomes of numerical calculations in a later section. The 
result makes sense because the mean population arising 
from a Poisson distribution with mean A at t = is equal 
to {4>)^ a= \. Intuitively, all initial situations sampled from 



such a distribution lead to a final population of unity as 
t — v oo except for the case of a initial population of zero, 
the probability of which is exp(— A). The probability of 
an initial situation that leads to a final population of 
unity is therefore 1 — exp(— A). The mean population 
as i — > oo is JV = x cxp(— A) + 1 x (1 — cxp(— A)) 
1 — exp(— A). 

We can confirm this result by considering kinetics 
starting from a definite initial population No- From Eq. 
( 35 1 we write 



V - #o^o JV °" 1 e MM (1 - exp(-^o)) 



n= 

^ SN n = 1 — <5jv o, 



71=1 



(70) 



n=l 



demonstrating that the asymptotic mean population is 
unity, unless Nq = in which case it is zero. 



I. Source enhanced coagulation 

We can determine the asymptotic behaviour of a sys- 
tem of coagulating species in which there is also an in- 
jection rate j. The underlying master equation is 



dP{N, t) 



jP(N-l,t)-jP(N,t)+ 



dt J ' K " J ' ' (71) 

k(N + 1)NP(N + 1, t) - kN(N - l)P(N, t), 

with the first term deleted for the case N = 0. The 
evolution of the corresponding pseudo-population is given 

by 



d(f> = jdt - K(f> 2 dt + i{2K) 1/2 cj)dWt. 

This time we find that 

d(e _< *) = -e~*d(t> - K^e^dt 
= -je'^dt - i(2n) 1/2 <i>e-' t 'dWt, 



(72) 



(73) 



so that d(e _ *) = -j(e _0 )dt and (e _</> ) oc e~ jt , with the 
proportionality factor depending on the initial condition. 
We also have 

d(<p n ) = nift 1 - 1 (j - K( f> 2 ))dt - K n(n - l)((/> n )dt, (74) 
which in a stationary state implies, for n > 1, that 



3 -i 

K 



in-l 



) st -(n-l)^") 6 



(75) 



This is reminiscent of a recurrence relatiorPS for modified 
Bessel functions, 7 n _i(x) — I n+ \{x) = (2n/x)I n (x), when 
written in the form %I n (x) — §I n —2(x) — (n— l)I n -i{x), 
which suggests that ((/>") st oc (j / K)"/ 2 /„_i(2(j / k) 1/2 ) for 



n > 0. Since with u = lwe have {<p" 
portionality factor is [/i(2(j'/k) 1 / 2 )] 



= j I ' k , the pro- 
Notice that we 
do not label these stationary averages according to c£>q 
since memory of the initial condition is lost for this case. 
In the final stationary state of a source enhanced coag- 
ulating system we therefore expect the mean population 
to be 



{4>h 



kJ h (2(j/«)V2) ! 



(76) 



which is similar to analysis performed elsewher e! 24 * 32 * 34 !. 
Note the limit of this result when j ' — > does not corre- 
spond to Eq. (|69| since dependence on the initial condi- 
tion is retained when j = 0. 



J. Remarks 

We conclude our consideration of analytical approaches 
by noting that the coagulation problem can be studied in 
a variety of ways, and that introducing complex stochas- 
tically evolving pseudo-populations, averaged over the 
noise and over a complex contour of initial values (or over 
the entire plane), might seem very elaborate compared 
with the direct analytic solution to the master equations, 
for example. Our purpose, however, is to establish explic- 
itly that the approach works for a simple case. Our atten- 
tion now turns to the numerical solution of the stochastic 
dynamics and averaging procedure to determine whether 
this can be performed efficiently for the model. Estab- 
lishing this would suggest that a wider set of stochastic 
problems based on master equations might be amenable 
to solution using these techniques. 



III. NUMERICAL APPROACH. 
A. Stochastic numerics 

An analytical solution to the SDEs corresponding to 
more complicated coagulation schemes, for example those 
involving multiple species and modelled in the mean field 
approximation with rate equations such as (|7| , is unlikely 
to be available. In these cases our approach is imple- 
mented through a numerical solution of the SDEs for the 
relevant pseudo-populations 0j(t), followed by averaging 
over the noise and a suitable weighting of the results over 
the initial values </>,[)• We test the feasibility of this strat- 
egy for the simple A + A ^> A model. 

We solve the SDE ( 15 1 numerically using a C++ code. 
In Figure [3j we plot Ke(4>)^, against t for (f) = 1, 2 
and 3, averaging 4> over 10 3 trajectories driven by in- 



dependent realisations of the Wiener process based on 
10 5 timesteps of length 5 x 10~ 5 , and with k = 1. Also 
shown are analytical results based on the formal series 
developed in section |IIE| truncated at the 12th power 




Figure 3. Numerical solutions to SDE ( 15 1, averaged over 10 3 
realisations of the noise with a timestep of 5 x 10 and k = 1, 
and with initial values <f>o = 1, 2 and 3, compared with series 
representations of (4>)<f> (dashed lines) including terms up to 
12th order in d>o, as derived in section III El 



of 4>o. This representation of the analytical solution ap- 
pears to be accurate for values of <j>$ up to three. For the 
largest value of t shown, the series approximates well to 
values 1 — exp(— (^o), in agreement with the analysis in 



section II H In contrast, he 12th order series representa- 
tion with 4>q = 5 deviates noticeably from the expected 
asymptotic behaviour for this range of t. indicating a se- 
ries truncation error. For small t the outcomes resemble 
the mean field result Eq. Q, as we suggest they should 
in Appendix A. 

However, the numerical results do not appear to be 
very satisfactory. Firstly, the averages are noisy, sug- 
gesting that 10 3 realisations are too few in number to 
obtain statistical accuracy, though this can of course be 
increased. Secondly, they do not seem to possess the cor- 
rect limits at large times: they all seem to tend towards 
unity. Thirdly, they often exhibit sharp peaks, which 
are sometimes large enough to cause the numerical sim- 
ulation to fail, as is seen in the case for <f>$ — 3, where 
a negative spike at t « 2.85 causes the calculation to 
crash. These instabilities, even if not terminal, have a 
disproportionate effect on the statistics when the num- 
ber of realisations is limited. Such insta bilities have been 
encountered before in simulation ;) 24 * 25 !, and in the next 
section we examine their origin and consider how they 
can be avoided. 



B. Elimination of instabilities 

The origin of the large fluctuations in the value of 4>(t) 
can be understood by considering the deterministic con- 
tribution to the evolution of <p(t) on the complex plane 
according to Eq. ( fl5"| . This is shown schematically by 
block arrows in Figure [4] representing the magnitude and 
orientation of the complex quantity — <f> 2 . Were it not for 
stochastic noise, a trajectory that started with a posi- 
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N, 58*- 



4 



4 



Re. 



ft 



k4 



Figure 4. The block arrows indicate the magnitude and direc- 
tion of drift of 4>{t) across the complex plane brought about 
by the deterministic term —K<j) 2 dt in the SDE (151. The ef- 



fect of this pattern is seen in an example realisation of the 
evolution. An excursion from the quasistable crescent shaped 
region near the origin towards negative Re(^) , and the return 
path in an anticlockwise direction, is responsible for a spike 
in the average of <f> shown in Figure [3] 



tive value of Re</> would remain in the right hand side 
of the complex plane and indeed would drift asymptot- 
ically towards the origin. The picture is quite different 
for a negative initial value of He(f>: there is a component 
of drift away from the origin, though in most cases the 
path ultimately moves into the right hand side of the 
plane and from there to the origin. The exception is 
for trajectories starting on the negative real axis, which 
asymptotically move towards — oo. This pattern of drift 
has the potential for instability. 

Of course stochastic dynamics is also driven by noise, 
and an example trajectory is illustrated in Figure [4] 
that for much of the time evolves in a rather well- 
behaved fashion over the complex plane, forming a cres- 
cent shaped trace situated to the right hand side of the 
origin. The noise maintains the trajectory away from 
the deterministic attractor at the origin. However, noise 
is also the undoing of this stability: if the trajectory 
wanders into the left hand side of the plane due to the 
stochastic term in Eq. (15), then a tendency to drift to- 
wards <f> = — oo can set in. This is also shown in Figure [1] 
at some point a fluctuation causes <j)(t) to move towards 
the left, away from the crescent pattern. In this example, 
quite a significant excursion towards negative real val- 
ues occurs before noise deflects the trajectory sufficiently 
away from the real axis, allowing the underlying drift to 
take <j>(t) in an anticlockwise fashion back towards the 
right hand side of the plane. Once there, the <j>(i) set- 
tles again into the previous pattern. The intervals spent 
within the crescent correspond to those periods in Figure 
[3] where fluctuations are small, whilst excursions into the 
left hand plane seem to be responsible for the spikes^. 



These instabilities are not desirable and action should 
be taken to eliminate them. They are less common if the 
timestep in the numerical simulation is reduced, since 
this lowers the likelihood of a noise-driven jump into the 
right hand side of the complex plane, but this increases 
the computational expense of the approach. We therefore 
need a mathematical scheme that can suppress the dan- 
gerous drift pattern in the stochastic dynamics while re- 
taining the statistical properties of solutions to the SDE. 

DrummoncP^ proposed a scheme for the elimination 
of such instabilities that took the form of a modification 
to the fundamental Poisson representation through the 
introduction of a weighting parameter f2(t), chosen to 
evolve according to 



dn = n g dw t , 



(77) 



with initial condition 0(0) = 1, and where g is an arbi- 
trary function. A new variable <p' is introduced, subject 
to the same initial condition imposed on <j>, and evolving 
according to 



d(j>' = -K(j)' 2 dt + i(2n) 1/2 (j)' (dW t - gdt) , 



(78) 



and Drummond showed that the f2-weighted average of 
(j)' over the noise is the same as the corresponding average 
of 0: 



(0) = w> 



(79) 



However, since <j)' and <f> evolve according to different 
SDEs, they need not suffer from the same instabilities. 

Drummond termed this the 'gauging away' of the in- 
stabilities of the original SDE, the terminology suggested 
by recognising that Eq. ( 79 1 possesses an invariance with 
respect to different choices of g. The approach is in fact 
equivalent to a transformation of the probability measure 
in stochastic calculus, and can be perhaps be most eas- 
ily understood as an application of the Cameron-Martin- 
Girsanov formulaPSlSSI j n stochastic calculus, as we de- 
scribe in more detail in Appendix B. 

We first explore one of Drummond's gauge functions 
g(4>') that has the capacity to tame the instability in the 
coagulation problem. Consider the choice 



g = i{K/2) 1/2 {tf-W\), 



which leads to the SDEs 



dtf = -K(f>'\<p'\dt + i(2K) 1/2 (f>'dW t . 



and 



dil = i(n/2) 



1/2 



) CldWt- 



(80) 



(81) 



(82) 



The crucial difference between Eqs. (811 and ( |15| is 
that the drift term in the SDE for <fi' is always directed 
towards the origin. It lacks the ability to produce an 
excursion towards <$>' = — oo. Meanwhile, £l(t) evolves 
diffusively in the complex plane, with no drift and hence 
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Re(Q(*)) 




Figure 5. Evolution of Re(0'S7) according to Eqs. (81 1 and 
Q82p, with cj>o = 1, 2 and 3, for the same timestep and num- 
ber of realisations as in Figure [3] The expected analytic be- 
haviour is indicated, and in comparison with Figure [3] the 
quality of the numerical results is much improved, which is 
a consequence of the gauging procedure. The evolution of 
Re(f2) for the case with cj}Q = 1 is shown in the inset, to illus- 
trate its slight deviation from the expected value of unity as 
t increases, due to sampling errors. 



d(Q) = 0. As long as <j>' is real, the increment dfi van- 
ishes. Because of this feature, Drummond called this 
choice of g a minimal gauge function, and regarded it as 
a natural choice for cases where <f>' is expected to take 
real values for most of its history. 

We have studied this reformulation of the coagulation 
problem and confirmed numerically that <fi'(t) does not 
suffer from instabilities of the kind experienced by <j)(t), 
and that the fi-weighted average of <f>' appears to agree 
with various analytic results expected for (<f>}$ , as shown 
in Figure [5j Nevertheless, we notice that the statistical 
uncertainty in the average of f2(t) grows as time pro- 
gresses, and this introduces a decline in accuracy, for a 
given number of realisations. We next address the rea- 
sons for this. 



C. Asymptotic behaviour of Q 

The diffusive behaviour of the weight function Q can 
be best demonstrated by considering the evolution of the 
mean of |fi| 2 . We have dfl = g$ldW t and hence dfi* = 
g*Vt*dWt such that, using Ito's lemma 



d\n\ 2 = ndn* + n*dn + g n g *n*dt 
= n g *n*dw t + n* g ndw t + \ g \ 2 \n\ 2 dt, 



and so 



d(\n\ 



(\ 9 \ 2 \n\ 2 )dt, 



(83) 



(84) 



which indicates that the mean square modulus of Q in- 
creases monotonically, whatever choice of function g is 




Reft 





early times 
late times 



Figure 6. Sketch of the probability distribution of Q over the 
complex plane as time proceeds. The common centre of the 
dark (early time) and lighter (late time) density plots lies at 
Q = 1. The diffusive behaviour suggests that the statistics 
in the sampling of Q and hence the quality of S7-weighted 
averages will suffer as time progresses. 



made. This result is illustrated schematically in Figure 
[6] the mean of Q is unity for all t, but its increasing 
mean square modulus suggests that the distribution of 
f2 spreads out. As a consequence, the extraction of the 
statistical properties of the problem will require an ever 
larger number of realisations as t increases. This is the 
price to pay, it seems, for the taming of the instabilities 
in the original SDE. 



D. Optimised gauge 

Nevertheless, the choices available to us in the gaug- 
ing procedure, or equivalently the freedom to shift the 
probability measure according to the Cameron-Martin- 
Girsanov formula, allow us to optimise the quality of the 
numerical estimates of (</>). Consider a new gauge func- 
tion labelled by the real variable R: 



g = t( K /2) 1/2 (0 



R), 



(85) 



which differs from Eq. (80 1 by an imaginary constant. It 
leads to the SDEs 



and 



dcj)' = -k4>' (R + \<t>'\)dt + ^nf^^'dWt 



<m = i (k/2) 1/2 ((/>'- 10' i - R) ndw t 



(86) 



(87) 



with the property (<f>) = (Q,<f>'). For R > this has the 
effect of strengthening the drift towards the origin in the 
SDE for <//. It is revealing to study the evolution of the 
square modulus of 4>' under this gauge. We write 



d(#Y*) = W* + + «(2k) 1/2 0' -i(2k) 1/2 0' 

= - K \(f>'\ 2 (R+\^\)dt~i(2 K ) 1 / 2 \(f ) '\ 2 dW t 

- k\c/)'\ 2 (R + \<f>'\) dt + i(2K) x l 2 \$\ 2 dW t + 2n\4>'\ 2 dt 



dt 



12 



= 2n{\- R)\cj)'\ 2 -2k\ 



A'|3 



(89) 



and thus | 0' | evolves deterministically towards an asymp- 
totic value of (1 — R), if R < 1, or zero if R > 1. 

It is also revealing to consider the SDE for f20' explic- 
itly: 

d(flcj)') = <j>'dD. + nd(j)' 

+ i (k/2) 1/2 (0' - |0'| - R) m^K^fidt 
= <t>'i (k/2) 1/2 (0' - |0'| - R) mW t 

+ n (-k4>' (R + \cf>'\) dt + i(2n) 1/2 (j)'dW t 

- K4>'n{<t>' - \4>'\ - R)dt 

= -K,n<f)' 2 dt + i (k/2) 1/2 0' (0' - |0'| - R + 2) QdW t , 

(90) 

such that as t — > oo the stochastic term tends towards 
either 

i (k/2) 1/2 0' ((1 - R) exp[i arg(0')] + 1) SldW t , (91) 

if R < 1, or i (k/2) 1/2 0' (2 - i?) QdW t if i? > 1. Further- 
more, the evolution of the average square modulus of f2 
is described by 
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Figure 7. Evolution of Re{Q(j>') gauged according to the im- 
proved procedure described by Eqs. ( 86 1 and ( 87 1 with R — 1, 
and for initial conditions fa = 1, 2 and 3. This should be 
compared with the results from Eqs. (pi} and 1 82 1 in Figure 



|5j obtained with the same timestep and number of realisa- 
tions but using R = 0. The evolution as t increases is less 
affected by statistical noise and corresponds more closely to 
the expected analytical behaviour. The evolution of f2 is more 
noisy, however, as indicated in the inset. 



d(\n\ 2 ) = -(\<f>> 



R\ 2 \n\ 2 )dt 



E. Further calculations 



(92) 



which tends to \n{\{\ - R) exp[i arg (<£')] - l\ 2 \n\ 2 )dt if 
R < 1, or \K,R 2 {\tt\ 2 )dt if R > 1. 

This analysis indicates that there are advantages in 
making specific choices of R. If we choose R = 2, Eq. 



Using the gauge function (85 1 with R = 1 we can 



( 90 1 shows that we can eliminate the noise term control- 



ling the asymptotic time evolution of Qcf)' . On the other 
hand, it is the diffusive behaviour of f2 that characterises 
the asymptotic quality of the statistics of other quanti- 
ties, so it is advisable to try to control this too. The 
strength of the diffusion depends on the quantity |(1 — 
i?)exp[iarg((//)]-l| 2 = (1 -R) 2 - 2(1 -R) cos arg(0') + 1 
for R < 1, and R 2 for R > 1. The behaviour is difficult 
to specify in detail, but for R > 1 it is clear that the 
asymptotic rate of diffusive spread of fi increases with 
R, such that large R is disadvantageous. 

We have explored gauges specified by values of R be- 
tween zero and two. As R increases, the diffusive be- 
haviour of becomes more marked while the noise in the 
evolution of £10' reduces. We find that the choice R = 1 
strikes a balance by suppressing noise in the evolution 
of Re(0)0 o = Ke(il(f>')^ , at the price of more significant 
statistical deviations of Re(il) from the expected value of 
unity. This is illustrated in Figure [7] which is to be com- 
pared with Figure [5j The shift in gauge seems to have 
transferred the statistical noise from (f2</>') into (f2). We 
suggest that R = 1 specifies an optimised gauge for this 
stochastic problem. 



study further aspects of the pseudo-population dynam- 
ics that were investigated analytically in section [TT] For 
particle coagulation starting from a Poisson distribution 
with mean A the averages are simply given by quantities 
A(N) = (A(4>))\. In Figure [i] we plot the evolution of 
a selection of state probabilities P(N,t) for A = 5, de- 
rived from Eq. ( 25 I , together with the solutions to Eq. 
([IJ for a similar case with the initial Poisson distribution 
truncated at N = 12 for convenience. The timestep and 
number of realisations are the same as in earlier cases. 
The correspondence is very good, and the noise does not 
significantly distort the results over the time scale for the 
completion of the coagulation. 

Finally, we examine source enhanced coagulation mod- 
elled by Eq. (72 1, with the use of the optimised gauge. 
The quantities (0) = (00') and (0 2 ) = (ft(/)' 2 ) are plotted 
for j = k = 1 in Figure [9] and they correspond well with 
the analytical results for late times obtained in section 
III II The initial condition is a Poisson distribution with 
A = 2. Furthermore, to indicate that not all quantities 
are statistically well behaved, we also show the average 
of exp(— 0), which ought to equal exp(— 2 — jt) for this 
case. Clearly more realisations would be needed to repro- 
duce this behaviour accurately with the prevailing choice 
of gauge. 
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t 



Figure 8. Probabilities of the occurrence of five, two and one 
particle in a coagulating system as time progresses, starting 
from a Poisson distribution with mean A = 5 at t = 0, com- 
paring numerical evaluations of Eq. ( |25[ ), gauged through Eq. 
(85 I with R = 1, with solutions to the master equations |T]). 
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Figure 9. Evolution of {cf>}^ and {(j> 2 }^ according to the 
source enhanced coagulation kinetics embodied in Eq. ( |72[ ), 
with j = k = 1, initial condition cj>o — A = 2, and 10 3 realisa- 
tions with time step 5 x 10 -5 . At late times they correspond 
well to the stationary state valu es (<ft ) Bt = /o(2)/ii(2) « 1.43 
b 2 ) s 



III 



The quantity (exp(- 



and {<f> } s t = 1 derived in section 

is also compared with the expected behaviour exp(— 2 
to illustrate that not every variable has acceptable statistics 
for this number of realisations of the stochastic evolution. 



it), 



IV. DISCUSSION 

We have outlined in this paper how to solve a simple 
coagulation problem using some rather elaborate analyt- 
ical and numerical methods, retrieving results that were 
previously either known, or computable from the master 
equations. The effort has not been misdirected, however, 
since our aim has been to establish that the approach 
can be implemented successfully for a well characterised 
example problem. Our real interest lies in more compli- 
cated cases of coagulation such as those with multiple 



species that agglomerate at a variety of rates. In such 
systems we expect to find behaviour that arises due to 
statistical fluctuations around low mean populations that 
cannot be captured in the usual mean field formulation. 
We plan to address such cases in further work. 

Specifically, we employ a Poisson representation^^ of 
the probability P(N, t) that there should be TV particles 
present at time t in the simple case of A + A — > A ki- 
netics. The Poisson representation is a superposition of 
Poisson distributions with complex means. The descrip- 
tion can be cast as a problem in the stochastic dynamics 
of a complex pseudo-population, the average of which 
over the noise and the initial condition is related to the 
average of the physical particle population. Analytical 
work gives a series expansion of this average in powers of 
the initial value of the pseudo-popu lation, and thi s can 
be employed to recover known resultP^ 28 30 32 34: . The 
development of the series involves the evaluation of mul- 
tiple integrals of functions of the Wiener process. We 
have also identified some exact results valid in the limit 
t — > oo, and in Appendix A give an approximate result 
for small nt that bears a resemblance to a mean field so- 
lution to the kinetics. We have discussed the evaluation 
of averages of general functions of the population, such 
as higher moments and the standard deviation. We have 
also evaluated moments of the pseudo-population in the 
stationary state of a coagulating system enhanced by a 
constant injection rate of new particles. 

Analytical work can only be performed in certain cir- 
cumstances, and more generally the numerical solution 
of the stochastic differential equations (SDEs) is neces- 
sary, followed by averaging over noise and initial con- 
dition. Unfortunately, numerical instabili ties ma ke this 
problematic, as has been noted previousl y 24 * 25 * 27 ! How- 
ever, we follow DrummoncP^in identifying an SDE for a 
pseudo-population that is free of instabilities, that is to 
be used together with a weighting procedure that repro- 
duces the statistics of the desired system. In Appendix 
B we have shown that this 'gauging' scheme is equiva- 
lent to an adaptive shift in the probability measure for 
the stochastic variable in the SDE, and that the weight- 
ing procedure is an application of the Cameron-Martin- 
Girsanov formula. This interpretation arguably makes 
the gauging procedure a little more intuitive. 

We have identified a choice of gauge, and associated 
transformation of the problem, that provides the statis- 
tical information in a rather optimised fashion. With 
this approach, the growth with time in the statistical 
uncertainty that is inherent with gauging can be ade- 
quately controlled, such that we need to generate rela- 
tively few realisations of the evolution, requiring a rather 
slight computational effort. A number of the analytically 
derived results have been reproduced using this numeri- 
cal approach. 

Of course Monte Carlo simulation of coagulation events 
taking place in an evolving population of particles would 
be an obvious alternative numerical method for studying 
this system. It is readily implemented, for example using 
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the Gillespie algorithirPH, and has the capacity to include 
population fluctuation effects, allowing us to explore var- 
ious processes of interest. It is an approach that is prob- 
ably easier to grasp than the methods we have outlined, 
and more computationally efficient for the problem under 
consideration here. Nevertheless, we believe that pseudo- 
population methods will prove to be more efficient than 
Monte Carlo in coagulation problems involving multiple 
species. 

Since Monte Carlo is equivalent to a solution of the 
master equations, the point can be made by enumerat- 
ing how many of these equations there might be. If we 
wish to study the evolution of clusters of monomers that 
have size dependent agglomeration properties, then each 
size must be treated as a distinct species. If we con- 
sider cluster sizes ranging from 1 to N max monomers, 
then the master equations will describe the evolution of 
probabilities P(Ni,N 2 ,--- >NN max ,t) where JVj denotes 
the population of clusters of size i. If the number of pos- 
sible values of each population is allowed to lie between 
and n™ ax , then the number of elementary population dis- 
tributions that must be considered is rij=i ax { n f lax + 
It is reasonable to say that this number could be rather 
large and will grow faster than iV max . The solution to 
such a large number of coupled ODEs, or the equivalent 
Monte Carlo simulation, would be daunting. 

On the other hand, the stochastic approach would 
require the solution to just N max SDEs for pseudo- 
populations 4>i(t); these are analogues of the mean pop- 
ulations of clusters of size i. The SDEs would bear some 
resemblance to the Smoluchowski coagulation equations 
for mean populations in the mean field, fluctuation-free 
limit, given in Eq. Q. Analytic solution to these equa- 
tions might not be possible, but numerical solution is not 
difficult, especially if we have techniques such as gauging 
at our disposal to avoid some of the pitfalls we have iden- 
tified. The drawback is that averaging of the SDEs over 
a variety of noise histories and perhaps initial conditions 
is necessary. Nevertheless, the task is linear in iV max and 
must eventually become more efficient than direct solu- 
tion to the master equations for large enough iV max . 

It is only when the mean of a population becomes small 
that deviations from mean field behaviour emerge, and 
so a hybrid approach might be possible whereby mean 
field rate equations are used to model the early stages of 
coagulation, going over to pseudo-population rate equa- 
tions when the mean population becomes small. This is 
possible since the mean field rate equations have a close 
resemblance to the evolution equations of the pseudo- 
populations, which is an important conceptual connec- 
tion. Specifically, we could integrate the equations for 
the <f>i(t) with the neglect of the noise when the modu- 
lus of 4>i{t) is large, starting them off at the initial mean 
populations of clusters of size i, such that they remain 
real, and introduce noise, and average over it, only when 
the modulus becomes small. This is to be investigated. 

In conclusion, we have made conceptual and numeri- 
cal developments of a method for kinetic modelling that 



was introduced some decades agcS3 but that appears not 
to have been fully exploited. We believe that in spite of 
some complexity in the formulation, the approach pos- 
sesses considerable intuitive value, and that it has the 
capability to treat systems with interesting statistical 
properties for which alternative methods are inappropri- 
ate or expensive. We intend to explore these possibilities 
in further studies. 
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Appendix A: Steepest descent evaluation 

We have seen that if we employ the representation Eq. 



( 29 1 , the set of initial values of O need only form a closed 



contour around the origin in the complex plane. We now 
demonstrate that for the purposes of computation there 
is an optimal contour, depending on the initial popula- 
tion in the problem N , and that 



(Al) 



where /at (</>o) is given by Eq. (Ill, reduces to Eq. Q 
for situations with a definite initial population, at least 
for small nt. 

Contour integrals of the form (fi exp[q(z)]dz may be 
evaluated by the method of steepest descents about 
saddle points in q[z) defined by locations z — z s in 
the complex plane where dq(z)/dz = 0. We repre- 
sent the integrand near a saddle point as exp(q(z)) 
exp(g(z s )) exp(q" (z s )(z — z s ) 2 /2) and integrate along a 
contour through z s on which the imaginary part of 
q"(z s )(z — z s ) 2 is zero and the real part is negative. If 
the contour is arranged such that it passes over saddle 
points along the line of steepest descent on the surface 
of the modulus of the integrand, and elsewhere follows 
valleys where the modulus of the integrand is small, then 
the contour integral may be approximated by a set of 
gaussian integrations about each saddle point position. 
We write 



N 



#0 7T^ ex P ' 
c Zm 



6 o -(A o + l)ln0 o +ln(0) 0o ), 



(A2) 

and consider nt -C 1 ■ In this limit M n ~ t n , according to 
Eq. (451, and hence Cj(t) s» (—nt)^ 1 such that ((f)) ^ ~ 
<f>o + 0(nt). Inserting this into the exponent in (A2|, 



we seek to identify saddle points in the remaining part, 
namely q((fj ) = (j> - N la^ . Solving q'(<p ) = 1 - 
No /4>o = identifies a single saddle point of the integrand 
at a location </>o = Nq. The modulus of 4>q N ° exp(</>o) has 
a minimum at this point and the real axis is the path 
of steepest ascent. The path of steepest descent passes 
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Re 0o 



exp((^o) for No = 5, indicating 



Figure 10. Modulus of 
the saddle point at (j>o = 5: the path of steepest descent 
through this point lies parallel to the imaginary axis. 



through the saddle point perpendicular to the real axis. 
The modulus of the integrand, for Nq = 5, is shown in 
Figure [10] 

Clearly, it is sensible to place the contour C along the 
path of steepest descent through the saddle point and 
then to complete the circuit around the origin through 
regions where the modulus of the integrand is as small as 
possible. The contour integral can then be approximated 
by the contribution along the straight line parallel to the 
imaginary axis through cf>Q ~ Nq. Writing c/)q = Nq + iy 
and with (4>)^, = </>o-F(</>o) we have 



N w exp(lnF(A )) 

i dy^ exp (q(N ) + ^"(AWo - N ) 2 )) 

Nq\, 



2tt 



(<f>)N N \ 



N °~\ 


— oo 


dye N °- N< - 


(Nq^ 


r 


exp(iVo) 






N N 0+ 1 



(0) 



A'n 



(A3) 

using Stirling's approximation n\ w n n cxp(— n)(27m) 1 / 2 
which is quite accurate even for Nq as low as unity. Fi- 
nally, the expansion ( 33 1 with the Cj (t) valid for nt <C 1 
may be written 



(A4) 



making Eq. (A3 1 and hence Eq. ( A2 1 consistent with the 
mean field approximation N w N /(l + nN t) for small 

Kt. 



Appendix B: Equivalence of Drummond gauging and 
the Cameron-Martin-Girsanov formula 



The equivalence between the statistical properties 
of the gauged stochastic variable <j>' , when suitably 
weighted, and _those of the ungauged vari able < />, derived 

is a con- 



by Drummond^ and explored in section III B 
sequence of some fundamental rules in stochastic calcu- 
lus that are expressed by the Cameron-Martin-Girsanov 
formulaPSIl An SDE such as 



dx = adt + bdWt . 



(Bl) 



is a statement of a connection between dx and a stochas- 
tic variable dWt with certain statistical properties. When 
we introduce expectation values such as (dx) we are im- 
plicitly defining a probability distribution or measure 
over the values taken by the variable dW t - Normally 
the notation dWt represents an increment in a Wiener 
process, the continuum limit of a symmetric 1-d random 
walk in discrete space and time, and the implication is 
that the probability distribution of values of dW t is gaus- 
sian with zero mean and variance equal to dt. 

But how might the expectation value of the increment 
dWt change if we were to evaluate it with respect to a 
different probability distribution? For example, what if 
it were distributed according to a shifted gaussian pro- 
portional to exp (— (dWt — m) 2 / (2dt)) where m is the 
non-zero mean of dWt making it no longer correspond to 
a Wiener process? Let us take the drift in the mean of 
dW t under such a shifted gaussian distribution to be pro- 
portional to the time elapsed, so we may write m = [idt 
and (dWtjQ — fidt. The new probability distribution, 
based on an asymmetric random walk and denoted Q, is 
indicated through a suffix on the expectation value. We 
write (dWtjp — in the old measure, which we denote 
P, and according to which dWt is indeed an increment 
in a Wiener process. 

Note that (dW t -^dt) Q = and ((dW t - ^dtf) Q = dt 
since the variance of dW t under measure Q is the same as 
that under P: we have shifted the gaussian probability 
distribution for dWt but not changed its width. Hence, 
if we define dW t = dW t — \xdt then (dWt)Q = and 
((dWt) 2 )Q = dt: we can identify a Wiener process that 
operates under probability measure Q, and relate it to a 
Wiener process under measure P. The SDE for x now 
reads 



dx = adt + b( dW t + lidt 



(B2) 



and we can choose to evaluate expectation values under 
measure P, for which (dWt)p = —fidt, or measure Q, for 
which (dW t ) Q = 0. 

The expectation value (dx)Q is the solution to a prob- 
lem that differs from the one initially posed, since the 
drift term in the SDE has been changed from adt to 
(a + bfi)dt. However, the point is that it is possible to es- 
tablish a link between expectation values under the two 
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different measures. The average of dx under measure P 
is equal to a weighted average of dx under Q. Formally, 
we can write 



(dx) 



dQ i q 



(B3) 



where dP/dQ is the Radon-Nikodym derivative of prob- 
ability measure P with respect to Q. Furthermore, the 
Cameron-Martin-Girsanov formula states that we can 
write 



dP 
dQ 



exp 



-ndW t 



1 



\?dt 



(B4) 



which is just a ratio of the gaussian distribu- 
tions exp[-(dW t + fxdt) 2 / (2dt)} and exp[-dW t 2 / (2dt)}. 
We now define a quantity Q through the relation 
exp(dmf2) = dP/dQ such that 



dlnQ = -[i 2 dt/2 - (j,dW t , 



(B5) 



and using Ito's lemma, we show that O evolves according 
to 

dil = d cxp (In Q) = exp (In CI) d In CI + - exp (In Cl) /j, 2 dt 
= Cl {-ii 2 dtj2 - fidWt + n 2 dt/2\ = -fifldW t . 

We see from Eqs. ( |B3p6| that averages of the incre- 
ment dx over differently distributed random increments 
can be related to one another, and that the quantity 
dCl describes the connection. In order for this to make 
sense mathematically, two conditions must be met. The 
first is that no process that has non-zero probability un- 
der P should be impossible under Q, and vice versa. 
The second, known as Novikov's condition, requires that 
(exp (i^ 2 di)) < oo. 

We now invert the point of view, and instead of con- 
sidering a single variable treated according to two dif- 
ferent averaging procedures, we relate the evolution of 
two different variables under the same averaging. Explic- 
itly, we consider variable x evolving as dx — adt + bdWt 
and another variable x' that evolves according to dx' = 
(a + bfi)dt + bdW t , with x(0) = x'(0). The above results 
imply that we can write 



(dx) = (exp(-fi 2 dt/2 - fidW t )dx') 



(B7) 



which resembles a combination of Eqs 
We shall demonstrate its validity short 



(|B3|} and <|B4| 



Note that 

there is no need to indicate the probability measure 
on the brackets, or further label the increment dWt, 
since here it is a standard Wiener increment with zero 
mean in the probability distribution under consideration. 
The inserted factor accounts for the difference in drift 
in the evolution of x and x and may also be written 
exp(dlnfi) with dlnfi = —fi 2 dt/2 — fidW t or equivalently 
dQ = -fiQdWf 



Clearly we can write exp(J dlnfi) = exp(lnO(i) 
lnf2(0)) = fi(f)/fi(0) so 

fi 2 (t')dt'- I ^(t')dWf ) . (B8) 



if we set Cl(0) = 1. Note that (Cl) = Cl(0) + J (dCl) = 
1 — J (fiCldWt) = 1. Novikov's condition on \l for a finite 

interval of time now reads (exp (i J* fi 2 (t')dt'^j 

Now consider the following: 



< oo. 



(Cl(t) (x'(t)-x'(0))) 

= (n e "^ 2(ti)dt_M(ti)dWi Z) d 4 



(B9) 



where we use a discrete representation of the time inte- 
gration. We have recognised that if k ^ j then 



/ e -kn 2 {tk)dt-^(t k )dw h 



TT^i^^dt-^dWi^ 



= / e -hn 2 (t k )dt-n(t k )dW k \ /TT e -y 2 (ti)dt-fi(ti)dWi dx ^ 



= /j] e -^ 2 ^) di -M*i)dm da .A 

Wfc / 

(BIO) 

such that by repetition of this step, the expression re- 
duces to (exp[— ^/j, 2 (tj)dt — n(tj)dWj]dx'j). 

Now, since (exp (— fidWt)) — exp([i 2 dt/2) we can write 

(exp {-fjidWt) dW t ) = -d (exp (-fJ,dW t )) /dfi 

= -fj,dt exp(^ 2 dt/2), (Bll) 

and with the insertion of dx' = (a + b^i)dt + bdWt, we 
conclude that 



s M 2 {tj )dt-n{tj )dW j dx r_ 



J l (B12) 
(a(tj) + bfi(tj)) dt — bfi(tj)dt = a(tj)dt = (dxj), 

such that 

(Cl(t)(x'(t) - x'(0))) = J2(dxj) = (a(t) - x(0)) , 

(B13) 

and since 

(Cl(t^'(0)) = (Cl(t))(x'(0)) = (x(0)), (B14) 
this means that 

(x(t)) = (Cl(t)x'(t)) . (B15) 

Thus we have shown that if we wish to evaluate the quan- 
tity (x(t)) generated by SDE dx = adt + bdWt, we could 
instead solve the SDEs 



dx' = adt + b(dW t - gdt) 
dQ = gCldW t , 



(B16) 
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with initial conditions x'(0) = x(0), 0(0) = 1, and with 
the function g(t) = —fJ,(t) taking arbitrary form subject 



to Novikov's condition, and then use the results to eval- 
uate the equivalent quantity (n(t)x'(t)). Clearly, this is 
identical to Drummond's gauging scheme. 
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